Cannot use getGEO because everything is in different places, use lumiBatch object to gather all the information
if(!file.exists('./../Data/Voineagu/LumiBatch_Voineagu.RData')){
create_Voineagus_pData_df = function(){
# Output: Dataframe with pData of Voineagu's raw data
# BrainBank, chip_array
Voineagu_suppMat = read.csv('../Data/Voineagu/aux/supplementary_table_2.csv')
Voineagu_suppMat$chip_array = paste(Voineagu_suppMat$Chip, Voineagu_suppMat$Array, sep='_')
Voineagu_suppMat = Voineagu_suppMat[,c('Brain.Bank.Case.Id','Cortex.area','Disease.status',
'chip_array')]
Voineagu_suppMat$Cortex.area = as.character(Voineagu_suppMat$Cortex.area)
# MATCH ALL FRONTAL + TEMPORAL IDS:
# GSM, BrainBank (relation obtained from the series matrix phenotype information)
ids_mapping = read.delim('../Data/Voineagu/aux/ids_mapping.txt', header=FALSE)
colnames(ids_mapping)[1:2] = c('GSM','Brain.Bank.Plus')
id_split = strsplit(as.character(ids_mapping$Brain.Bank.Plus),'_')
ids_mapping$Disease.status = sapply(id_split, `[`, 1)
ids_mapping$Brain.Bank.Case.Id = sapply(id_split, `[`, 2)
ids_mapping$Cortex.area = sapply(id_split, `[`, 3)
cortex_area = list('C'='cerebellum','T'='temporal','F'='frontal')
ids_mapping$Cortex.area[ids_mapping$Cortex.area %in% names(cortex_area)] =
cortex_area[ids_mapping$Cortex.area[ids_mapping$Cortex.area %in% names(cortex_area)]]
disease_status = list('A'='autism', 'C'='control')
ids_mapping$Disease.status[ids_mapping$Disease.status %in% names(disease_status)] =
disease_status[ids_mapping$Disease.status[ids_mapping$Disease.status %in% names(disease_status)]]
# BrainBank, metadata
samples_metadata = read.delim2('../Data/Voineagu/aux/supplementary_table_1.tsv')
# Frontal and temporal info:
ids_mapping_f_t = merge(Voineagu_suppMat, ids_mapping,
by=c('Brain.Bank.Case.Id','Disease.status','Cortex.area'))
# MATCH ALL CEREBELLUM IDS:
# GSM, BrainBank, chip_array (cerebellum)
ids_mapping_c = read.delim('../Data/Voineagu/aux/cerebellum_ids_mapping.txt', header=FALSE)
colnames(ids_mapping_c) = c('GSM','Brain.Bank.Plus','chip_array')
id_split = strsplit(as.character(ids_mapping_c$Brain.Bank.Plus),'_')
ids_mapping_c$Disease.status = sapply(id_split, `[`, 1)
ids_mapping_c$Brain.Bank.Case.Id = sapply(id_split, `[`, 2)
ids_mapping_c$Cortex.area = sapply(id_split, `[`, 3)
ids_mapping_c$Cortex.area[ids_mapping_c$Cortex.area %in% names(cortex_area)] =
cortex_area[ids_mapping_c$Cortex.area[ids_mapping_c$Cortex.area %in% names(cortex_area)]]
ids_mapping_c$Disease.status[ids_mapping_c$Disease.status %in% names(disease_status)] =
disease_status[ids_mapping_c$Disease.status[ids_mapping_c$Disease.status %in% names(disease_status)]]
# JOIN FRONTAL+TEMPORAL INFO WITH CEREBELLUM AND ADD METADATA:
ids_mapping_all = rbind(ids_mapping_f_t, ids_mapping_c)
ids_mapping_all = ids_mapping_all[,c('GSM','chip_array','Brain.Bank.Case.Id','Cortex.area',
'Disease.status')]
samples_full_data = merge(ids_mapping_all, samples_metadata, by='Brain.Bank.Case.Id')
samples_full_data$Cortex.area = as.character(samples_full_data$Cortex.area)
# Add age group data
age_group = cut(as.numeric(samples_full_data$AGE), c(-0.5, 0, 0.6, 10, 20, 30, 40, 50, 60, 70, 80),
labels=c('Fetal','Infant','Child','10-20','20s','30s','40s','50s','60s','70s'))
samples_full_data$age_group = as.character(age_group)
return(samples_full_data)
}
# AssayData data
# Supplementary file GSE28521_non-normalized_data.txt.gz
LumiBatch_v = lumiR('../Data/Voineagu/GSE28521_non-normalized_data.txt.gz')
# Phenotype data
phenoData = create_Voineagus_pData_df()
phenoData = phenoData[match(colnames(exprs(LumiBatch_v)), phenoData$chip_array),]
phenoData$SEX[phenoData$SEX=='5'] = 'F'
rownames(phenoData) = phenoData$GSM
# Change subject IDs from chip_array to GSM format
chipArray_to_gsm = phenoData$GSM
names(chipArray_to_gsm) = phenoData$chip_array
#colnames(LumiBatch_v) = chipArray_to_gsm[colnames(LumiBatch_v)]
# LumiBatch object
pData(LumiBatch_v) = phenoData
save(LumiBatch_v, file='../Data/Voineagu/LumiBatch_Voineagu.RData')
remove(assayData, phenoData, chipArray_to_gsm)
} else {
load(paste('../Data/Voineagu/LumiBatch_Voineagu.RData', sep='/'))
datExpr = exprs(LumiBatch_v)
datMeta = pData(LumiBatch_v)
}
# Get Biomart information
ensembl = useMart("ENSEMBL_MART_ENSEMBL",dataset="hsapiens_gene_ensembl", host="feb2014.archive.ensembl.org")
getinfo = c("illumina_humanref_8_v3", "ensembl_gene_id","external_gene_id", "entrezgene",
"chromosome_name", "start_position", "end_position")
geneDat = getBM(attributes = getinfo, filters="illumina_humanref_8_v3",
values = rownames(datExpr), mart = ensembl)
idx = match(rownames(datExpr), geneDat[,"illumina_humanref_8_v3"])
datGenes = cbind(rownames(datExpr), geneDat[idx,])
rownames(datGenes) = datGenes[,1]
datGenes$length = datGenes$end_position-datGenes$start_position
# SFARI Genes
SFARI_genes = read_csv('./../Data/SFARI/SFARI_genes_with_ensembl_IDs.csv')
# GO Annotations
GO_annotations = read.csv('./../Data/GO_annotations/genes_GO_annotations.csv')
GO_neuronal = GO_annotations %>% filter(grepl('neuron', go_term)) %>%
mutate('ID' = as.character(ensembl_gene_id)) %>%
dplyr::select(-ensembl_gene_id) %>% distinct(ID) %>%
mutate('Neuronal' = 1)
# Combine SFARI and GO information
gene_info = data.frame('ID'=datGenes$ensembl_gene_id) %>% left_join(SFARI_genes, by='ID') %>%
mutate(`gene-score`=ifelse(is.na(`gene-score`), 'None', `gene-score`)) %>%
left_join(GO_neuronal, by='ID') %>%
mutate('gene.score'=ifelse(`gene-score`=='None' & Neuronal==1, 'Neuronal', `gene-score`)) %>%
mutate('gene.score'=ifelse(is.na(gene.score), 'None', gene.score))
rm(ensembl, geneDat, GO_annotations, LumiBatch_v, getinfo, idx, create_Voineagus_pData_df)
Remove cerebellum samples:
datMeta = datMeta[datMeta$Cortex.area!='cerebellum',]
datExpr = datExpr[,datMeta$chip_array]
RNA-Seq for 58 cortical brain-tissue samples across frontal and temporal lobes, comprising 29 samples from control subjects and 29 samples from ASD subjects
print(paste0('Dataset includes ', nrow(datExpr), ' genes from ', ncol(datExpr), ' samples belonging to ', length(unique(datMeta$Brain.Bank.Case.Id)), ' different subjects.'))
## [1] "Dataset includes 24526 genes from 58 samples belonging to 32 different subjects."
Sex distribution: There are five times more males than females
table(datMeta$SEX)
##
## F M
## 9 49
Age distribution: Subjects between 5 and 56 years old with a mean close to 29
summary(datMeta$AGE)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 5.00 18.25 28.50 28.79 39.00 56.00
Diagnosis distribution: Fairly balanced
table(datMeta$Disease.status)
##
## autism control
## 29 29
Brain region distribution:
table(datMeta$Cortex.area)
##
## frontal temporal
## 32 26
Diagnosis and brain region are very well balanced
table(datMeta$Disease.status, datMeta$Cortex.area)
##
## frontal temporal
## autism 16 13
## control 16 13
1. Filter genes with start or end position missing
to_keep = !is.na(datGenes$length)
datGenes = datGenes[to_keep,]
datExpr = datExpr[to_keep,]
#rownames(datGenes) = datGenes$ensembl_gene_id
print(paste0('Removed ', sum(!to_keep), ' genes, ', sum(to_keep), ' remaining'))
## [1] "Removed 780 genes, 23746 remaining"
2. Filter genes with low expression levels
The density distribution doesn’t seem to suggest any clearly defined threshold. Selecting 300
plot_data = data.frame('id'=rownames(datExpr), 'mean_expression' = rowMeans(datExpr))
ggplotly(plot_data %>% ggplot(aes(x=mean_expression)) +
geom_density(color='#0099cc', fill='#0099cc', alpha=0.3) +
geom_vline(xintercept=300, color='gray') + scale_x_log10() +
ggtitle('gene Mean Expression distribution') + theme_minimal())
to_keep = rowMeans(datExpr)>300
datGenes = datGenes[to_keep,]
datExpr = datExpr[to_keep,]
print(paste0('Removed ', sum(!to_keep), ' genes, ', sum(to_keep), ' remaining'))
## [1] "Removed 8943 genes, 14803 remaining"
Using node connectivity as a distance measure, normalising it and filtering out genes farther away than 2 standard deviations from the left (lower connectivity than average, not higher)
Using \(s_{ij}=|bw(i,j)|\) to define connectivity between genes.
Filtering 4 autism samples
absadj = datExpr %>% bicor %>% abs
## alpha: 1.000000
netsummary = fundamentalNetworkConcepts(absadj)
ku = netsummary$Connectivity
z.ku = (ku-mean(ku))/sqrt(var(ku))
plot_data = data.frame('ID'=rownames(datMeta), 'sample'=1:length(z.ku), 'distance'=z.ku, 'Region'=datMeta$Cortex.area, 'Diagnosis'=datMeta$Disease.status)
ggplotly(plot_data %>% ggplot(aes(sample, distance, color=Diagnosis)) + geom_point() + geom_hline(yintercept=-2, color='gray') + theme_minimal())
print(paste0('Outlier samples: ', paste(as.character(plot_data$ID[plot_data$distance< -2]), collapse=', ')))
## [1] "Outlier samples: GSM706448, GSM706425, GSM706447, GSM706455"
to_keep = abs(z.ku)<2
datMeta = datMeta[to_keep,]
datExpr = datExpr[,to_keep]
print(paste0('Removed ', sum(!to_keep), ' samples, ', sum(to_keep), ' remaining'))
## [1] "Removed 4 samples, 54 remaining"
rm(absadj, netsummary, ku, z.ku, plot_data, to_keep)
print(paste0('After filtering, the dataset consists of ', nrow(datExpr), ' genes and ', ncol(datExpr), ' samples'))
## [1] "After filtering, the dataset consists of 14803 genes and 54 samples"
# Transform Signal to Ratio
#row_mean = apply(datExpr, 1, mean)
#ratio = (datExpr)/(row_mean+1)
#exprs(LumiBatch_v) = as.matrix(ratio[,colnames(ratio) %in% colnames(datExpr)])
#exprs = as.matrix(ratio[,colnames(ratio) %in% colnames(datExpr)])
#LumiBatch_v = ExpressionSet(assayData = exprs)
LumiBatch_v = ExpressionSet(assayData = datExpr)
rownames(datMeta) = datMeta$chip_array
pData(LumiBatch_v) = datMeta
# Variance Stabilisation
LumiBatch_v = lumiT(LumiBatch_v, method = 'log2', ifPlot = TRUE)
## Perform log2 transformation ...
# Normalisation
LumiBatch_v = lumiN(LumiBatch_v, method='rsn')
## Perform rsn normalization ...
## 2019-09-03 23:13:45 , processing array 1
## 2019-09-03 23:13:45 , processing array 2
## 2019-09-03 23:13:46 , processing array 3
## 2019-09-03 23:13:46 , processing array 4
## 2019-09-03 23:13:46 , processing array 5
## 2019-09-03 23:13:46 , processing array 6
## 2019-09-03 23:13:46 , processing array 7
## 2019-09-03 23:13:46 , processing array 8
## 2019-09-03 23:13:47 , processing array 9
## 2019-09-03 23:13:47 , processing array 10
## 2019-09-03 23:13:47 , processing array 11
## 2019-09-03 23:13:47 , processing array 12
## 2019-09-03 23:13:47 , processing array 13
## 2019-09-03 23:13:47 , processing array 14
## 2019-09-03 23:13:47 , processing array 15
## 2019-09-03 23:13:47 , processing array 16
## 2019-09-03 23:13:47 , processing array 17
## 2019-09-03 23:13:47 , processing array 18
## 2019-09-03 23:13:47 , processing array 19
## 2019-09-03 23:13:47 , processing array 20
## 2019-09-03 23:13:47 , processing array 21
## 2019-09-03 23:13:47 , processing array 22
## 2019-09-03 23:13:47 , processing array 23
## 2019-09-03 23:13:48 , processing array 24
## 2019-09-03 23:13:48 , processing array 25
## 2019-09-03 23:13:48 , processing array 26
## 2019-09-03 23:13:48 , processing array 27
## 2019-09-03 23:13:48 , processing array 28
## 2019-09-03 23:13:48 , processing array 29
## 2019-09-03 23:13:48 , processing array 30
## 2019-09-03 23:13:48 , processing array 31
## 2019-09-03 23:13:48 , processing array 32
## 2019-09-03 23:13:48 , processing array 33
## 2019-09-03 23:13:48 , processing array 34
## 2019-09-03 23:13:48 , processing array 35
## 2019-09-03 23:13:48 , processing array 36
## 2019-09-03 23:13:48 , processing array 37
## 2019-09-03 23:13:49 , processing array 38
## 2019-09-03 23:13:49 , processing array 39
## 2019-09-03 23:13:49 , processing array 40
## 2019-09-03 23:13:49 , processing array 41
## 2019-09-03 23:13:49 , processing array 42
## 2019-09-03 23:13:49 , processing array 43
## 2019-09-03 23:13:49 , processing array 44
## 2019-09-03 23:13:49 , processing array 45
## 2019-09-03 23:13:49 , processing array 46
## 2019-09-03 23:13:49 , processing array 47
## 2019-09-03 23:13:49 , processing array 48
## 2019-09-03 23:13:49 , processing array 49
## 2019-09-03 23:13:49 , processing array 50
## 2019-09-03 23:13:49 , processing array 51
## 2019-09-03 23:13:49 , processing array 52
## 2019-09-03 23:13:50 , processing array 53
## 2019-09-03 23:13:50 , processing array 54
save(LumiBatch_v, file=paste0('../Data/Voineagu/LumiBatch_v_preprocessed.RData'))
# Update datExpr rownames
rownames(datExpr) = datGenes[rownames(datExpr),]$ensembl_gene_id
datExpr = datExpr[!is.na(rownames(datExpr)),]
datExpr = log2(datExpr)
mod = model.matrix(~Disease.status, data=datMeta)
corfit = duplicateCorrelation(datExpr, mod, block=datMeta$Brain.Bank.Case.Id)
fit = eBayes(lmFit(datExpr, mod,block=datMeta$Brain.Bank.Case.Id,
correlation = corfit$consensus),trend=T)
DE = topTable(fit,coef=2, number=Inf, sort.by="none")
plot_data = DE %>% left_join(gene_info, by='ID')
ggplotly(plot_data[abs(plot_data$logFC)<1,] %>% ggplot(aes(x=gene.score, y=abs(logFC), fill=gene.score)) +
geom_boxplot() + scale_fill_manual(values=SFARI_colour_hue(r=c(1:6,8,7))) +
theme_minimal() + theme(legend.position='none') + xlab('Gene score') + ylab('abs(lfc)'))
No clear relation between SFARI score and mean expression/SD. This dataset doesn’t seem to have the same problem as Gandal’s
plot_data = data.frame('ID'=rownames(datExpr), 'MeanExpr'=rowMeans(datExpr),
'SDExpr'=apply(datExpr,1,sd)) %>%
left_join(gene_info, by='ID')
p1 = ggplotly(plot_data %>% ggplot(aes(gene.score, MeanExpr, fill=gene.score)) + geom_boxplot() +
scale_fill_manual(values=SFARI_colour_hue(r=c(1:6,8,7))) + theme_minimal() +
theme(legend.position='none') + xlab('Gene score') + ylab('mean expression'))
p2 = ggplotly(plot_data %>% ggplot(aes(gene.score, SDExpr, fill=gene.score)) + geom_boxplot() +
scale_fill_manual(values=SFARI_colour_hue(r=c(1:6,8,7))) + theme_minimal() +
ggtitle('Mean Expression (left) and SD (right) by SFARI score') +
theme(legend.position='none') + xlab('Gene score') + ylab('standard deviation'))
subplot(p1, p2, nrows=1)
rm(plot_data, p1, p2)
Data seems to be close to homoscedastic
plot_data = data.frame('ID'=rownames(datExpr), 'Mean'=rowMeans(datExpr), 'SD'=apply(datExpr,1,sd))
plot_data %>% ggplot(aes(Mean, SD)) + geom_point(color='#0099cc', alpha=0.1) +
geom_smooth(method='lm', color='gray', size=0.5) + theme_minimal()
rm(plot_data)
No Batch information available in the metadata!
PCA: There doesn’t seem to be a specific pattern related to diagnosis or brain region
*There weren’t any patterns in MDS and t-SNE either (MDS was the same as PCA)
pca = datExpr %>% t %>% prcomp
datMeta$ID = datMeta$chip_array
plot_data = data.frame('ID'=colnames(datExpr), 'PC1' = pca$x[,1], 'PC2' = pca$x[,2]) %>%
left_join(datMeta, by='ID') %>%
dplyr::select('PC1','PC2','AGE','age_group','SEX','Disease.status', 'Cortex.area')
p1 = plot_data %>% ggplot(aes(PC1, PC2, color=Disease.status)) + geom_point() + theme_minimal() +
ggtitle('PCA') +
xlab(paste0('PC1 (',round(100*summary(pca)$importance[2,1],1),'%)')) +
ylab(paste0('PC2 (',round(100*summary(pca)$importance[2,2],1),'%)'))
p2 = plot_data %>% ggplot(aes(PC1, PC2, color=Cortex.area)) + geom_point() + theme_minimal() + ggtitle('PCA') +
xlab(paste0('PC1 (',round(100*summary(pca)$importance[2,1],1),'%)')) +
ylab(paste0('PC2 (',round(100*summary(pca)$importance[2,2],1),'%)'))
grid.arrange(p1, p2, nrow=1)
rm(pca, plot_data)
First Principal Component explains 96.6% of the total variance
There’s a really strong correlation between the mean expression of a gene and the 1st principal component
pca = datExpr %>% prcomp
plot_data = data.frame( 'PC1' = pca$x[,1], 'PC2' = pca$x[,2], 'MeanExpr'=rowMeans(datExpr))
plot_data %>% ggplot(aes(PC1, PC2, color=MeanExpr)) + geom_point(alpha=0.5) + theme_minimal() +
scale_color_viridis() + ggtitle('PCA') +
xlab(paste0('PC1 (',round(100*summary(pca)$importance[2,1],1),'%)')) +
ylab(paste0('PC2 (',round(100*summary(pca)$importance[2,2],1),'%)'))
rm(pca, plot_data)
save(datExpr, datMeta, datGenes, file='./../Data/Voineagu/preprocessed_data.RData')
#load('./../Data/Voineagu/preprocessed_data.RData')
sessionInfo()
## R version 3.6.1 (2019-07-05)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 18.04.3 LTS
##
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/openblas/libblas.so.3
## LAPACK: /usr/lib/x86_64-linux-gnu/libopenblasp-r0.2.20.so
##
## locale:
## [1] LC_CTYPE=en_GB.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_GB.UTF-8 LC_COLLATE=en_GB.UTF-8
## [5] LC_MONETARY=en_GB.UTF-8 LC_MESSAGES=en_GB.UTF-8
## [7] LC_PAPER=en_GB.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_GB.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] parallel stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] viridis_0.5.1 viridisLite_0.3.0 gridExtra_2.3
## [4] plotly_4.9.0 forcats_0.4.0 stringr_1.4.0
## [7] dplyr_0.8.2 purrr_0.3.2 readr_1.3.1
## [10] tidyr_0.8.3 tibble_2.1.3 ggplot2_3.2.0
## [13] tidyverse_1.2.1 WGCNA_1.68 fastcluster_1.1.25
## [16] dynamicTreeCut_1.63-1 limma_3.40.6 biomaRt_2.40.4
## [19] lumi_2.36.0 GEOquery_2.52.0 Biobase_2.44.0
## [22] BiocGenerics_0.30.0
##
## loaded via a namespace (and not attached):
## [1] tidyselect_0.2.5 robust_0.4-18.1
## [3] RSQLite_2.1.2 AnnotationDbi_1.46.1
## [5] htmlwidgets_1.3 grid_3.6.1
## [7] BiocParallel_1.18.1 munsell_0.5.0
## [9] codetools_0.2-16 preprocessCore_1.46.0
## [11] statmod_1.4.32 nleqslv_3.3.2
## [13] withr_2.1.2 colorspace_1.4-1
## [15] knitr_1.23 rstudioapi_0.10
## [17] stats4_3.6.1 robustbase_0.93-5
## [19] labeling_0.3 GenomeInfoDbData_1.2.1
## [21] bit64_0.9-7 rhdf5_2.28.0
## [23] vctrs_0.2.0 generics_0.0.2
## [25] xfun_0.8 R6_2.4.0
## [27] doParallel_1.0.15 GenomeInfoDb_1.20.0
## [29] illuminaio_0.26.0 locfit_1.5-9.1
## [31] bitops_1.0-6 reshape_0.8.8
## [33] DelayedArray_0.10.0 assertthat_0.2.1
## [35] promises_1.0.1 scales_1.0.0
## [37] nnet_7.3-12 gtable_0.3.0
## [39] affy_1.62.0 methylumi_2.30.0
## [41] rlang_0.4.0 zeallot_0.1.0
## [43] genefilter_1.66.0 splines_3.6.1
## [45] rtracklayer_1.44.3 lazyeval_0.2.2
## [47] acepack_1.4.1 impute_1.58.0
## [49] broom_0.5.2 checkmate_1.9.4
## [51] BiocManager_1.30.4 yaml_2.2.0
## [53] modelr_0.1.5 crosstalk_1.0.0
## [55] GenomicFeatures_1.36.4 backports_1.1.4
## [57] httpuv_1.5.1 Hmisc_4.2-0
## [59] tools_3.6.1 nor1mix_1.3-0
## [61] affyio_1.54.0 RColorBrewer_1.1-2
## [63] siggenes_1.58.0 Rcpp_1.0.1
## [65] plyr_1.8.4 base64enc_0.1-3
## [67] progress_1.2.2 zlibbioc_1.30.0
## [69] RCurl_1.95-4.12 prettyunits_1.0.2
## [71] rpart_4.1-15 openssl_1.4
## [73] bumphunter_1.26.0 S4Vectors_0.22.0
## [75] SummarizedExperiment_1.14.1 haven_2.1.1
## [77] cluster_2.1.0 magrittr_1.5
## [79] data.table_1.12.2 mvtnorm_1.0-11
## [81] matrixStats_0.54.0 mime_0.7
## [83] hms_0.5.1 evaluate_0.14
## [85] xtable_1.8-4 XML_3.98-1.20
## [87] mclust_5.4.5 readxl_1.3.1
## [89] IRanges_2.18.2 compiler_3.6.1
## [91] minfi_1.30.0 KernSmooth_2.23-15
## [93] crayon_1.3.4 htmltools_0.3.6
## [95] later_0.8.0 mgcv_1.8-28
## [97] pcaPP_1.9-73 Formula_1.2-3
## [99] rrcov_1.4-7 lubridate_1.7.4
## [101] DBI_1.0.0 MASS_7.3-51.4
## [103] Matrix_1.2-17 cli_1.1.0
## [105] quadprog_1.5-7 GenomicRanges_1.36.0
## [107] pkgconfig_2.0.2 fit.models_0.5-14
## [109] GenomicAlignments_1.20.1 registry_0.5-1
## [111] foreign_0.8-72 xml2_1.2.2
## [113] foreach_1.4.7 annotate_1.62.0
## [115] rngtools_1.4 pkgmaker_0.27
## [117] multtest_2.40.0 beanplot_1.2
## [119] XVector_0.24.0 bibtex_0.4.2
## [121] rvest_0.3.4 doRNG_1.7.1
## [123] scrime_1.3.5 digest_0.6.19
## [125] Biostrings_2.52.0 rmarkdown_1.13
## [127] base64_2.0 cellranger_1.1.0
## [129] htmlTable_1.13.1 DelayedMatrixStats_1.6.0
## [131] curl_3.3 shiny_1.3.2
## [133] Rsamtools_2.0.0 nlme_3.1-141
## [135] jsonlite_1.6 Rhdf5lib_1.6.0
## [137] askpass_1.1 pillar_1.4.2
## [139] lattice_0.20-38 httr_1.4.0
## [141] DEoptimR_1.0-8 survival_2.44-1.1
## [143] GO.db_3.8.2 glue_1.3.1
## [145] iterators_1.0.12 bit_1.1-14
## [147] stringi_1.4.3 HDF5Array_1.12.2
## [149] blob_1.2.0 latticeExtra_0.6-28
## [151] memoise_1.1.0